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Reionization by UV or X-ray sources 
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Abstract. We present simulations of the 21-cm signal during the epoch of reionization. We focus on properly 
modeling the absorption regime in the presence of inhomogeneous Wouthuysen-Field effect and X-ray heating. We 
ran radiative transfer simulations for three bands in the source spectrum (Lyman, UV, and X-ray) to fully account 
for these processes. We find that the brightness temperature fluctuation of the 21 cm signal has an amplitude 
greater than 100 mK during the early reionization, up to 10 times greater than the typical amplitude of a few 10 
mK obtained during the later emission phase. More importantly, we find that even a rather high contribution from 
QSO-like sources only damps the absorption regime without erasing it. Heating the IGM with X-ray takes time. 
Our results show that observations of the early reionization will probably benefit from a higher signal-to-noise 
value than during later stages. After analyzing the statistical properties of the signal (power spectrum and PDF) 
we find three diagnostics to constrain the level of X-ray, hence the nature of the first sources. 



1. Introduction 

The epoch of reionization (EoR) started with the forma- 
tion of the first sources of light around z = 15 — 30. As 
show n by the Gunn-Peterson effect ( Gunn fe Peterson! 
Il965h in the spectra of high-Fredshift quasars (QSO) 
(e.g., iFan et all 12006). the universe was fully reionized 
by z ~ 6. WMAP 5-year results show that the optical 
depth for the Thomson scattering of CMB photons trav- 
eling through the reio nizing universe is r = 0.084 ± 0.016 
(jKomatsu et al l l2009h . Together with the Gunn-Peterson 
data, this strongly favors an extended reionization period 
between z > 11 and z = 6. 

While other observa tions, such a s the L yman-a emit- 
ter luminosity function ( Ouchi et al.l . 20091) . may produce 
other constraints on the history of reionization in the 
next few years, the most promising is the observation 
of the 21-cm line in the neutral IGM using large radio- 
interferometers (LOFAB0, MWA@, GMRTQ 21-CMA0, 
SKA0). The signal will be observed in emission or in 
absorption against the CMB continuum. Both theoreti- 



cal modeling (Madau ct al. ■ H997l:lFurlanetto et all 



and simulation s (e.g-JCiardi fe M adau 2003; Gnedin 



Mellema et"aH l2006bl: iLidz et al.l 120081 llchikawa et al 



2006) 



2004; 



2009; iThomas et al.l 120091 ) show that the brightness tem- 
perature fluctuations of the 21 cm signal have an ampli- 
tude of a few 10 mK in emission, on scales from tens of 
arcmin down to sub-arcmin. With this amplitude, and 
ignoring the issue of foreground cleaning residuals, sta- 
tistical quantities such as the three-dimensional power 
spectrum should be measurable with LOFAR or MWA 
with a few 100 hours integration ([Morales Sz Hewitt , 2004 ; 
Furlanetto et all . 120061: ILidz et a!T~2008 ). In absorption 
however, the amplitude of the fluctuations may exceed 
100 mK (iGnedinl . 12004 ISantos et all l2008t iBaek et al 
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2009), the exact level depending on the relative contribu- 
tion of the X-ray and UV sources to the process of cosmic 
reionization. 

The signal will be seen in absorption during the initial 
phase of reionization, probably at z > 10, when the accu- 
mulated amount of emitted X-ray is not yet sufficient to 
raise the IGM temperature above the CMB temperature. 
The duration and intensity of this absorption phase, reg- 
ulated by the spectral energy distribution (SED) of the 
sources, are crucial. SKA precursors able to probe the rel- 
evant frequency range, 70 - 140 MHz, may benefit from a 
much higher signal-to-noise than during later periods in 
the EoR. However, if the absorption phase is confined at 
redshifts above 15, RFI and the ionosphere will become 
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an problem. Quite clearly, the different types of sources 
of reionization and different formation histories produce 
very different properties for the 21-cm signal. It is there- 
fore important for future observations to explore the range 
of astrophysically plausible scenarios using numerical sim- 
ulations. 

To properly model the signal, it is necessar y to use > 
100/i _1 Mpc box sizes (jBarkana fc Loebl . |2003) . Together 
with a large box size, it is desirable to resolve halos with 
masses down to 10 8 M Q as these contain sources (able to 
cool below their virial temperature by atomic processes), 
or even minihalos with masses down to 10 4 Mq became 
they act as an efficie nt photon s ink be cause of their high 
recombination rate ( Iliev et all . l200otl . As this work fo- 
cuses on improving the physical modeling, we restrict our- 
selves to resolving halos with a mass 10 10 Af Q or higher. 
Indeed, simulating the absorption phase correctly, as we 
do in this work, requires a more extensive and more costly 
implementation of radiative transfer. We are exploring the 
direct implication of this improved physical modeling, and 
will turn to better mass resolution in the near future. 

There are three bands in the sources SED that in- 
fluence the level of the 21 cm signal: the Lyman band, 
the ionizing UV band, and the soft X-ray band. Lyman 
band photons are necessary to decouple the spin temper- 
ature of hydrogen from the CMB temperature through 



the W outhuysen-Field effect (|Wouthuvsenl . Il952t [Field 



19581 ). and make the EoR signal visible. UV band pho 



tons arc of course responsible for the ionization of the 
IGM, and soft X-rays are able to preheat the neutral gas 
ahead of the ionizing front, deciding whether the decou- 
pled spin temperature is less (weak preheating) or greater 
(strong preheating) than the CMB temperature. While 
a proper modeling should perform the full 3D radiative 
transfer in all 3 bands, a simpler modeling has often been 
used in previous works. Indeed, for the usual source SEDs 
and source formation histories, once the average ioniza- 
tion fraction of the universe is greater than ~ 10%, the 
background flux of Lyman-a photons is so high that the 
hydrogen spin temperature is fully coupled t o the kinetic 



temperature by the Wouthuysen-Field effect (Baek et al 



l2009f h Thereafter, computing the Lyman band radiative 
transfer is unnecessary. In the same spirit, it has usually 
been assumed that the preheating of the IGM by soft 
X-ray was strong enough to raise the kinetic tempera- 
ture much higher than the CMB temperature everywhere 
early in the EoR. However, both assumptions fail during 
the early reionization: the absorption phase. Even in the 
later part of reionization the second assumption may fail, 
depending on the nature of the sources. We will quan- 
tify this possibility in this paper. Computing the full ra- 
diative transfer in all three bands is necessary to study 
the absorption regime. Indeed, fluctuations in the local 
Lyman-a flux induce fluctuations in the spin tempera- 
ture (while the Wouthuysen-Field effect is not yet sat- 
urated), which, i n turn, modify the power spectrum of 
the 2 1 cm signal (Barkana fc Loebl 120051: ISemelin et al. , 
l2007t IChuzhov fc Zhentd . l2007t iNaoz fc Barkanai 12008: 



Baeket al.l [20091) . The same is true for the fluctuations in 
the local flux of X - ray p hotons (jPritchard fc Furlanettd . 



20071: ISantos et all 120081) 



Let us emphasize however that, in modeling Lyman- 



a and X-ray fl uctua ti ons, [B arkana fc Loe 



Naoz fc Barkanai (l2008h IPritchard fc Furlanettd (|200 



and ISantos et al.l (|2008h all use the semi-analytical 
approximation that the IGM has a uniform density 
of absorbers and ignore win g effects in the rad iative 



transfer of Lyman - a pho tons. ISemelin et al.l (|2007l ) and 



Chuzhov fc Zheng (|2007l) have shown that these wing ef- 



fects do exist. Moreover, once reionization is under way, 
ionized bubbles create sharp fluctuations in the number 
density of absorbers (not to mention simple matter density 
fluctuations). In this work, for the first time, we present 
results based on simulations with full radiative transfer 
for both Lyman-a and X-ray photons. 

What are the possible candidates as sources of reion- 
ization? Usually, two categories are considered: ionizing 
UV sources (Pop II and III stars), and X-ray sources 
(quasars). When we study 21 cm absorption, however, we 
must distinguish between Pop II and Pop III stars be- 
yond the large difference in luminosity per unit mass of 
formed star. Indeed Pop II stars have a three times larger 
Lyman band to ionizing UV band luminosity ratio than 
Pop III stars. It means that the 21 cm signal will reach 
its full power (near saturated Wouthuysen-Field effect) at 
a lower average ionization fraction for Pop II stars than 
for Pop III stars. The relevant question is: how long do 
Pop III stars dominate the source population before Pop 
II stars take over? The answer to this question, related to 
the whole process of star formation, feedback and metal 
enrichment of the IGM, is a difficult one. At this stage, 
state of the art numerical simulations of the EoR use sim- 



ple prescriptions in the best case (e.g. llliev et al.ll2007fl . or 
simply ignore this issue. 

The other category of sources are X-ray sources. They 
may be mini-quasars. X-ray binaries, supernovae (jOhl , 
20011: iGlover fc Brand) . 2003), or even more exotic can 



didates such as dark stars (jSchleicher et al. . 20091 ). The 
exact level of emission from these sources is a matter 
of speculation. The generally accepted view is that stars 
dominate ov er X-ray sources and are sufficient to drive 
reion i zation (|Shapiro fc Girouxl 119871: iGiroux fc Shapiro! 
I l996l iMadau et all 119991 ICiardi et all l2003h . Recently. 
IVolonteri fc Gnedinl (|2009l) supported the opposite view. 
While, in their models, X-ray sources are marginally able 
to complete reionization by z ~ 6, they fi nd a very low 



contri bution from stars. Indeed they rely on lGnedin et al 



(2008J) who find, using numerical simulations, a negligi- 
ble escape fraction for ionizing radiations from galaxies 
with total mass less than a few 1O 1O M , who should actu- 
ally contribute t o 90% of the ionizing photon production 
during the EoR ( Choudhurv fc Ferraral 120071) . While the 
physical modeling in their innovative simulations is quite 
detailed, this surprizing behavior of the escape fraction 
definitely needs to be checked at higher resolution and 
with different codes. For the time being the best Simula- 
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tions can only explore a plausible range of X-ray contri- 
butions, and quantify the impact on observables. When 
the observations become available we would like to be 
able, using simulation results, to derive tight constraints 
on the relative level of emission from ionizing UV and X- 
ray sources. This work, exploring the 21 cm signal for a 
few different levels of X-ray emission, is a first step toward 
this goal. 

The paper is organized as follows. We present the nu- 
merical methods in §2 and describe our source models in 
§3. In §4, we show the results and analyze the differences 
between the models. We discuss our findings and conclude 
in §5. 

2. Numerical simulation 

The numerical me thods used i n this work are similar to 
those presented in iBaek et al. | (|2009h (hereafter Paper I). 
The references to previous and some new validation tests 
are presented in the Appendix. The dynamical sim ula- 
tions have been run with GADGET-2 (jSpringell [20051) and 
post-processed with UV continuum radiative transfer and 
further processed with Ly-a transfer using LICORICE. 
The same cosmological parameters and particle number 
are used and we refer the reader to Paper I for details re- 
lated to the numerical methods and parameters. The main 
improvement on the previous work is using a more real- 
istic source model including soft X-ray and implementing 
He chemistry. 

We have run seven different simulations, all of which 
use the same lOOft^Mpc box, density fields, and star 
formation rate, but with different initial mass functions 
(IMF), chemistry (with helium or without), X-ray frac- 
tion of the total luminosity or X-ray spectral index. SI 
is the reference model. S2 has a top-heavy IMF (Salpeter 
IMF restricted to a 100 - 120M© range), while the oth- 
ers uses a Salpeter IMF in a 1.6 — 120Af Q range. Only S3 
contains helium. In all other models, helium is replaced 
by the same mass of hydrogen. X-ray radiative transfer is 
included in S4, S5, S6 and S7. They have either different 
X-ray fraction of the total luminosity or X-ray spectral 
index. The basic parameters of these simulations are sum- 
marized in Tab. [1] 

The simulations are controlled by a few parameters. 
We adopted the same value as in paper I for the maxi- 
mum value of the number of particles per radiative trans- 
fer cell in the adaptive grid: N max — 30. The resulting 
minimum radiative transfer cell size is 200 hr x kpc at 
z = 6.6. Between two snapshots, i. e. ~ 10 Myr, we cast 
3 x 10 6 photon packets for photoionization (all in the UV 
for models SI to S3, half in the UV and half X-rays for 
models S4 to S7), and 3 x 10 7 photons for Lyman-a trans- 
fer. At the end of the simulations (z ~ 6), the number of 
sources reaches ~ 15000, so the number of ionizing photon 
packets per source is only 200. However, at this final stage 
the sources are highly clustered and very large and ion- 
ized regions surround the source clusters. So the clustered 
sources cooperate to reduce the Monte Carlo noise at the 



ionization fronts. In addition, the adaptive grid re- 
sponds better than a fixed grid to sampling issues: 
big cells where there are few photons, sm al l cells 
where there are many. iMaselli &: Ferrara (2005) 
presents convergence tests for a Monte-Carlo ra- 
diative transfer code very similar to ours. Their 
convergence tests suggest that the typical level of 
noise in our ionization and temperature cubes is 
~ 10%. We accept is as a reasonable value, espe- 
cially since, having run the Ilie v et al.l ( 20091 ) com- 
parison tests, we are confident that our ionization 
fronts propagate at the correct speed. We use 1000 
frequency bins in each of the photoionizing-UV and X-ray 
spectra. For Lyman-a transfer, we sample the frequency 
at random between Lyman-a and Lyman- f3. 



Model 




IMF 


Helium 


-kstar 


Lqso 


spectral index 


SI 


1.6 


- 120M Q 


No 


100 % 


% 




S2 


100 


- 120M© 


No 


100 % 


% 




S3 


1.6 


- 120M© 


Yes 


100 % 


% 




S4 


1.6 


- 12OM 


No 


99.9 % 


0.1 % 


a=1.6 


S5 


1.6 


- 120M© 


No 


99.9 % 


0.1 % 


«=0.6 


S6 


1.6 


- 12OM 


No 


99 % 


1 % 


a=1.6 


S7 


1.6 


- 120M Q 


No 


90 % 


10 % 


a=1.6 



Table 1. Simulation parameters. L star is the stellar lumi- 
nosity fraction and Lqso is the X-ray luminosity fraction 
of the total luminosity. 



2.1. X-ray radiative transfer 

The main difference between the cosmological radiative 
transfer of ionizing UV and X-ray is the mean free path 
of the photons, at most a few 10 comoving Mpc in the 
first case, possibly several 100 Mpc in the second case. A 
usual trick in implementing UV transfer is to use an infi- 
nite speed of light: do so with LICORICE (see Paper I). 
This is correct if the crossing time of the photons between 
emission and absorption points is much less than the re- 



comb ination time, the photo-ionization time (jAbel et al 
1999h and the typical time for the variation of luminos- 
ity of the sources. This is the case in most of the IGM 
during the EoR, except very close to the sources where 
the photo-ionizing rate is very high. Obviously, this is not 
the case for X-rays which have a much longer crossing 
time. Consequently, we implemented the correct propaga- 
tion speed for X-ray photons. We propagate an X-ray pho- 
ton packet during one radiative transfer time step At reg 
(< 1 Myr, same notation from Fig.l of Paper I) over a dis- 
tance of cAt reg , where c is the velocity of light. Then, the 
frequency of the photon packet and photoionizing cross 
sections are recomputed with the updated value of the 
cosmological expansion factor. The photon packet prop- 
agates during the next radiative transfer time step using 
these new parameters. If the photon packet loses 99% of its 
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initial energy, we drop it. X-ray photon packets containing 
photons with an energy of several keV pass through the pe- 
riodic simulation box several times before they lose most of 
their initial content. For each density snapshot, that is ev- 
ery ~ 10 Myr, 1.5 millions of photon packets are sent from 
the X-ray sources. About half of them are absorbed during 
the computation on the same density snapshot when they 
were emitted, and the other half is stored in memory to 
be propagated through the next density snapshots. Some 
X-ray packets with very high energy photons still survive 
several snapshots later, so the number of stored photon 
packets grows as simulations progress. About 50 millions 
photon packets are stored in memory toward the end of 
the simulations. 

It may seem that this memory overhead, which sets 
a limitation to the possible simulations with LICORICE, 
would not appear with radiative transfer algorithms which 
naturally include a finite velocity of light like moment 
methods. However, these methods would suffer from an 
overhead connected to the number of frequency bins nec- 
essary to correctly model X-rays, while it does not exist in 
Monte-Carlo methods. Including complete X-ray transfer 
in EoR simulations comes at a non- negligible cost, what- 
ever the numerical implementation. Since the X-ray pho- 
tons can propagate over several box sizes during several 
tens of Myr, the X-ray frequency can redshift consider- 
ably between emission and absorption. The cross-section 
of photoionization has a strong frequency dependence, 
so we have to redshift the frequency of the photons. At 
each radiative transfer time step At reg , we update the 
frequency of all the X-ray photon packets, 



V{t + At reg ) 



a(t) 



a(t + At reg ) 



(1) 



3. Source model 

3.1. Computing the star formation rate 

Our new source model needs the star formation rate for 
all baryon particles. We recompute the star formation in 
the radiative transfer simulations rather than to rerun the 
dynamical simulation. Here is why and how. 

We adopted the procedure described in 



Mihos fc Hernauistl ( 19941 ). employing a local Schmidt 
law and an hybrid-particles algorithm to implement it 
in our code. Indeed, in our model, the star formation 
rate solely depends on the density, and we make the 
assumption that the star formation feedback (kinetic and 
thermal) on the dynamics does not vary much from the 
fiducial simulation. LICORICE uses the classical Schmidt 
law: 

dU 1 



-j- = — (if Pg > ^threshold), 



where f » is defined by : 



t* — tn 



Pg 



^threshold 



-1/2 



(2) 



(3) 



p g is the gas density and /* is the star fraction. We set the 
parameters to* and ^threshold so that the evolution of the 
global star fraction follows closely that of the S20 simula- 
tion (20ft. -1 Mpc) in Paper I, and reionization completes 
at z ~ 6. In this way, we reuse the tuning made for the S20 
simulation, and at high z, we get a similar star formation 
history as in the higher resolution (but smaller box size) 
S20 simulation. All simulations in the present work have 
a 100/i -1 Mpc box size. Following the above equations, 
a gas particle whose local density exceeds the threshold 
(^threshold = 5 ^critical x H&) increases its star fraction, /*, 
where p C riticai is the critical density of the universe and fib 
is the cosmological baryon density parameter. 



where a(t) is the expansion factor of the Universe. 

The treatment of non-thermal electrons produced by 
X-ray will be described in ^ 13.51 

2.2. Helium reionization 

The intergalactic medium is mainly composed of hydro- 
gen and helium, with contributions of 90% and 10% in 
number. Until now, we have run simulations with hydro- 
gen only, but including helium is worth studying because 
the different value of the ionization thresholds and pho- 
toionization rates could affect the reionization history. We 
included He, He +, and He + + in LICORICE, and used 
ai (Il992l) and IVerner et al. (|l996t) for various cooling 
rate and cross sections. When helium chemistry is turned 
on, the ionization fractions (H + , He + and Hc ++ ) and the 
temperature are integrated explicitly using the adaptive 
scheme described in Paper I. More details on the numer- 
ical methods and a validation tests of the treatment of 
helium arc presented in Appendix. 



3.2. Limiting the number of sources 

We compute the increase in the star fraction for each par- 
ticle, A/* between two consecutive snapshots. Then the 
total mass of young stars formed in a particle is m x A/* , 
where m is the mass of the particle. To avoid a huge 
number of sources, we had to set a threshold on the new 
star fraction for the particle to act as a source. We used 
A/* > 0.001. We checked that this leaves out a negligible 
amount of star formation, about 0.4%. It happens that 
several source particles reside in the same radiative trans- 
fer cell, but we treated individually ray tracing for each 
source. 



3.3. Choosing an IMF 

With our mass resolution, this amount m gas x A/» corre- 
sponds to a star cluster or a dwarf galaxy so an IMF should 
be taken into account. We choose a Salpeter IMF, with 
masses in the range 1.6M Q - 120M Q or 100M© - 120M Q 
(model S2). The first range is used to model the SED of 
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an intermediated Pop II and Pop III star population, and 
the other one is for pure Pop III stars. 

3.4. Computing the luminosity and SED of the stellar 
sources 

The next step is to make the link between the amount 
of created stars and the luminosity of the sources. When 
only the ionizing UV luminosity is considered, it is quite 
justified to use simple models. For example we can make 
it proportional to the mass of the host dark matter halo 
( Iliev et all 2006bl) . or, as we did in Paper I, to the mass of 
the baryonic particles newly converted into stars. Things 
are more complicated when we consider both the Lyman 
band and ionizing UV. Indeed, since each particle is mas- 
sive enough to contain a representative sample of the 
choosen IMF, and since each mass bin has a different 
life time, we should consider an SED evolving with the 
age of the star particle. This would be possible using pre- 
tabulated SEDs. However, unlike in Paper I, we decided 
to use hybrid particles which begin to produce photons as 
soon as a small fraction of the particle is turned into stars. 
This is useful to make the local luminosity less noisy in 
the early EoR when the source mass resolution is an issue. 
Including this star formation history for each particle and 
convolving with the time- varying SED would be extremely 
costly in terms of both memory and computation time. 

We simplified the issue by considering the fact that 
in the Lyman and UV band, most of the luminosity is 
produced by the massive stars, with a short life time com- 
parable to the time between two snapshots of our sim- 
ulations. So we decided to use a constant SED and lu- 
minosity during a characteristic life time. Both luminos- 
ity and SED are computed independently in each of the 
Lyman and ionizing band. To compute the luminosity and 
the SED for a star particle we use the data for massive, 
low metallicity (Z = 0.004 Zp,) stars in the ma i n sequ ence 



mass [M e ] log(L/L M@ ) log(T cff ) tn te [Myr] 



( Mevnet fe Maederl . 12003: Hansen fe Kawaleit 119941 ) (see 
Tab. [2]). The details of how this is done can be found in 
Appendix. The constant luminosity and characteristic life 
time, computed in the two spectral bands, are given in 
Tab. [3] We find characteristic life time of < 8 Myr for 
the UV band. In the implementation of the UV transfer 
however, for technical reasons, the source fraction of the 
particle actually shines for a duration equal the interval 
between two snapshots. This varies varies between 6 and 
20 Myr, so we recalibrate the luminosity to produce the 
correct amount of energy. The whole point of the proce- 
dure, is to take the different typical life time in the Lyman 
band into account, especially at at z, when Lyman-a cou- 
pling is not yet saturated. We should not concentrate the 
emission within a single snapshot interval, which is 3 times 
shorter than the source life-time, or we would artificially 
boost the coupling between t he spin temperature of hy- 
drogen and the kinetic temperature of the gas and alter 
the resulting brightness temperature. Consequently we let 
each newly formed star fraction of a particle shine for 3 



120 


6.3 


4.7 


3 


60 


5.8 


4.6 


4.5 


40 


5.6 


4.5 


6 


30 


5.2 


4.5 


7 


20 


4.8 


4.45 


10 


15 


4.65 


4.4 


14 


12 


4.2 


4.37 


20 


9 


3.8 


4.3 


34 


5.9 


2.92 


4.18 


120 


2.9 


1.73 


3.97 


700 


1.6 


0.81 


3.85 


3000 



Table 2. Physical properties of low metalicity (Z = 
0.004Z Q ) main sequence stars ( Mevnet fe Maeder . 20051) 
L is the bolometric luminosity, (T c g) is the effective tem- 
perature and tnfe [Myr] is the life time of the star. 



consecutive snapshots, which is close to the typical life 
time in the Lyman band, and we still recalibrate the lu- 
minosity to produce the correct amount of energy. While 
we do not use a time-evolving SED, we believe that imple- 
menting different life times for the Lyman and UV sources 
with the correctly average luminosities is a substantial im- 
provement in our source model. 

We use an escape fraction f esc = 0.12 for photoionizing 
UV photons and f esc = 1 for Lyman-a. 




Fig. 1. Normalized spectral intensity of our source model. 
Black solid line is the SED from Salpeter IMF and red 
dashed line is from top-heavy IMF. 



3.5. X-ray source model 

X-rays can have a significant effect on the 21 cm brightness 
temperature. The X-ray photons, having a smaller ioniz- 
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IMF mass range 


Energy band 


10.24 eV < E < 13.6 eV 


E > 13.6 eV 


1-120 [M Q ] 


Luminosity [erg/s]^ 


6.32 x 10 44 


2.14 x 10 45 




Life time[Myr] A 


20.36 


8.03 


100-120 [M Q ] 


Luminosity [erg/s] s 


9.96 x 10 45 


3.12 x 10 46 




Life time[Myr] s 


3.32 


3.31 



Table 3. Averaged luminosities and life times of our source model for a baryon particle depending on the energy 
band. A values are from Salpeter IMF and B values are form top-heavy IMF. 



ing cross-section, can penetrate neutral hydrogen further 
than UV photons and heat the gas above the CMB tem- 
perature. This X-ray heating effect on the IGM is often 
assumed to be homogeneous because of X-rays' long mean 
free path. In reality, the X-ray flux is stronger around 
the sources and the inhomogeneous X-ray flux can bring 
on ex tra fluctuations for the 21 cm brightness tempera - 



ture (jPritchard fc Furlanettol . 120071: ISantos et all 120081 ). 



Moreover, patchy reionization induce further fluctuations 
in the local X-ray flux which can only be accounted for 
using a full radiative transfer modeling. 

X-rays luminosity 

First, we need to determine the luminosity and location 
of X-rays sources. We simply divided the total luminos- 
ity, Ltot-, °f all source particles into a stellar contribution, 
L st ar, and a QSO contribution, Lqso- Lqso depends on 
the star formation rate, since L to t itself is proportional 
to the increment of the star fraction, A/*, between two 
snapshots of the dynamic simulation. 

One reason for this approach is that X-ray binaries and 
supernova remnants contribute to X-ray sources as well as 
quasars and they are strongly r elated to the star forma - 
tion rate. Following the work of Gl over fc Brand! (l2008h . 
we took 0.1% of L t ot as the fiducial X-ray source lumi- 
nosity, Lqso- However, considering that they assumed a 
simple and empirically motivated model we have also run 
simulations with different values Lqso, 1% and 10% of 
Ltot- Quasar luminosity fractions less than 0.1% are not 
of interest for us, since their heating effect will be negligi- 
ble. 



X-ray energy range and nature of the sources 

First, we have to choose the photon energy range since 
hard X-ray photons have a huge mean free path which 
costs a lot for ray-tracing computations. Th e comoving 
mean free path of an X-ray with energy E is ( Furlanettol 
200fih 



X x = 4.9x~l /3 



l + z 



( E 



, . comoving Mpc. 

15 J \300eV J 8 y 

(4) 

Only photons with energy below E ~ 2[(1 + 
z)/15] 1 /%H / I 3 fceV are absorbed within a Hubble time and 
the E~ 3 dependence of the cross-section means that heat- 



ing is dominated by soft X-rays, which do fluctuate on 
small scales (Furlanetto et al. 2006). Therefore, we choose 
an energy range for X-ray photons from O.lkeV to 2fceV. 
The photons with energy higher than 2/ceV are not ab- 
sorbed until the end of simulation at z ~ 6. 

While the most likely astrophysical sources of X-ray 
during the EoR are supernovae, X-binaries and (mini-) 
quasars, it is interesting to mention that the X-ray SED 
of su pernovae and X-binaries typically peaks above 1 fceV 
(e.g. IOhll200lh . This means that most of the X-rays emit- 
ted by these sources will interact with the IGM more 
than 10 s years later, which is not true for QSO-like SEDs. 
During this time interval the global source mass (and, to 
first order, luminosity) easily rises by a factor of 10. Thus 
the longer delay will lower the effective luminosity of X- 
binaries and supernovae compared to QSO. For this rea- 
son, but also to avoid detailed modeling of some aspects 
while others, like the overall luminosity of X-ray sources, 
remain largely unconstrained, we use QSO as our typical 
X-ray source. 

QSO spectral Index 

We model the specific luminosity of our QSO-like sources 
as a power-law with index a; 



L v = k\ — 



A; is a normalization constant so that 

r v 2 



rV2 

Lqso = L v Av 



(5) 



(6) 



where hv\ = O.lfceV and hv\ = 2keV. The amount of 
X-ray heating can be altered by the shape of the spec- 
trum but there exists a large observational uncertainty in 
the mean and distribution of a. We extrapol ate from the 



measu rement of extreme UV spectral index by Tel fer et al 



(2002) to higher energy. The index values are derived from 
fitting lRy < E < iRy, and we extrap olated to 2fceV. The 
measured value bv lTelfer et al. ( 20021) is approximately ~ 
1.6, but with a lar ge gaussian standard deviation of 0.86. 
IScott et af . ( 20041 ) derived an average value of a = 0.6 
from a sample of FUSE and HST quasars. We choose 
a = 1.6 as our fiducial case, and used a = 0.6 for compar- 
ison. 
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Secondary Ionization 

X-rays deposit energy in the IGM by photoionization 
through three channels. The primary high velocity elec- 
tron torn from hydrogen and helium atoms distributes its 
energy by 1) collisional ionization, producing secondary 
electrons, 2) collisional excitation of H and He and 3) 
Coulom b collision with therma l electr ons. The fitting for- 
mula in Shull fc van Steenbera (1985) is used to compute 
the fraction of secondary ionization and heating. These 
are taken into account when computing the evolution of 
the state of the IGM. 

Then, it is legitimate to ask whether the lyman-a elec- 
trons resulting from the collisional excitations are impor- 
tant for the Wouthuysen-Field effect. The simple answer is 
that, in our choice of models, the energy emitted as X-ray 
is at most 10% of the UV energy, itself 3 times less than the 
Lyman luminosity. Moreover at m ost 40% of the X-ray en- 
ergy is converted into excitations ( Shull fe van Steenberd . 
1985). and only ~ 3 0% of the excitations r esult in a 
Lyman-alpha photon ( Pritchard fc Furlanetto . 2006). So, 
in the best case, Lyman-a photons produced by X-ray rep- 
resent only 0.3% of the photons produced directly by the 



4. Results 

4.1. Ionization fraction 

The evolution of the averaged ionization fraction tells us 
about the global history of reionization. We plot the mass 
and volume weighted average ionization fraction in Fig. [2] 
Including quasar (S4, S5, S6 and S7, not plotted) does 
not change the global evolution of the ionization fraction 
much, because of its small fraction of the total luminos- 
ity. The total number of emitted photons is similar to SI. 
The S3 simulation has also the same number of emitted 
photons, but S3 reaches the end of reionization a little bit 
earlier (Az « 0.25) than SI. Unlike other simulations, S3 
contains helium which occupies 25% of the IGM in mass. 
Including helium, the total number of atoms is reduced 
by 20%. At the same time, unless X-ray are the dominant 
source of reionization, most of the helium is only ionized 
once while z > 6, due to the higher energy threshold for 
the secondary ionization(Hc ++ ). Therefore the number of 
emitted photons per baryon is higher in S3 than SI, and 
it results in an earlier reionization. 

On the other hand, S2 has the same number of photon 
absorbers as SI, but the total number of emitted photons 
is much higher than in SI. Using a top-heavy IMF, it 
produces 10 times more photons (see Fig. [1] and Tab. EJ), 
and results in a Az « 1 earlier reionization. 

In all three cases, volume weighted values are less than 
mass weighted values, since gas particles in dense regions 
around the sources are ionized first. The volume occupied 
by each particle is estimated using the SPH smoothing 
length. 

We computed the Thomson optical depth for all sim- 
ulations, the values are r = 0.062, 0.076, 0.064 for SI, S2 



and S3. The other simulations (S4-S5) have the same r 
as SI, since they follow the same evolution of ionization 
fraction. They are somewhat lower t han the Thompson oi 



tical depth derived from WMAP5 (jHinshaw et all . 1200 



t = 0.084 ± 0.016, only the S2 value is within la of the 
WMAP5 value. A variable escape fraction, decreasing with 
time, would allow the IGM to start ionization earlier and 
increase t, without terminating ionization after z=6. 



S2 : Mass average 

52 : Volume average 

53 : Mass average 
S3 : Volume average 
SI : Mass average 
SI : Volume average 




Fig. 2. Mass (thicker lines) and volume (thinner lines) 
weighted ionization fraction of hydrogen. SI and S3 use 
Salpeter IMF and S2 use top-heavy IMF. S3 contains he- 
lium clement. 



4.2. Gas temperature 

The main goal of this study is to investigate the effect of 
inhomogeneous X-ray heating on the 21-cm signal. If the 
Ly-a coupling is sufficient, and X-rays can heat the gas 
above the CMB temperature Tcmb, the 21-cm signal will 
be observed in emission. However if the X-ray heating is 
not very effective, particularly during the early phase of 
the EoR, we will observe the signal in absorption. 

We plot in Fig. [3] the averaged gas temperature of the 
neutral IGM whose ionization fraction xhii is less than 
0.01. We chose the criterion of xhii < 0.01 for the fol- 
lowing reasons. Once a gas particle is 10% ionized, it is 
heated by photoheating to a temperature of several thou- 
sand Kelvin. At redshift 10, the number of gas particles 
which have an ionization fraction between 0.01 and 0.1 
is only 0.1% of all the particles, but if we include these 
particles, the average temperature increases from 2.94K 
to 5.41K. Therefore, we used the criterion xhii < 0.01 
to evaluate properly the average temperature of neutral 
regions, and verified that xhii < 0.001 gives a very sim- 
ilar average temperature. We have checked that even for 
model S7 which has the highest level of X-rays, at z > 7.5, 
99% of the neutral IGM has indeed an ionization fraction 
less than 1%, so we have not excluded a significant frac- 
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tion of the 21 cm emitting IGM from our average. This 
neutral gas is mostly located in the voids of the IGM. In 
fact, we have to consider Ly-a heating as well as X-ray 
heating since a few K difference can reduce the intensity 
of 8Tb by up to 100 mK. We recompute the gas temper- 
ature to include Ly-ct heating as a post-treatm ent using 
the formula from Furlanetto fc Pritchardl (<2006f l . This was 
detailed in Paper I. The temperature of all simulations in 
Fig. [3] decreases until z w 12 because of the adiabatic 
expansion of the universe. Then S7, which has the high- 
est Lqsoi starts to increase first and reaches the CMB 
temperature at redshift z w 8.8. Our fiducial model, S4, 
which contains 0.1% of total energy as X-rays, shows very 
little increase with respect to SI, a simulation without 
X-rays. Even for S7, the gas temperature of neutral hy- 
drogen in the void is still under the Tqmb until z « 8.8. 
This means that the X-ray heating needs time to heat the 
IGM above Icmb- Even with a rather high level of X- 
ray, the absorption phase survives and produces greater 
brightness temperature fluctuations than the subsequent 
emission phase (the delay in the absorption of X-ray con- 
nected to the long mean free path is partly responsible 
for this). It will be important to keep this result in mind 
when choosing the design and observation strategies for 
the future instruments. 



CMB 

S7 : 10 % quasar 
S6 : 1% quasar 
S5 : 0.1% quasar 
S4: 0.1% quasar 
SI : no quasar 




10 

Redshifl z 



Fig. 3. The evolution of the gas temperature of the neutral 
IGM with redshift. The neutral gas is chosen so that its 
ionization fraction is less than 0.01. 



There is a large observational uncertainty in the mean 
and distribution of the quasar spectral index a. Our fidu- 
cial model assumes a = 1.6 but we run a simulation (S5) 
with a = 0.6 for comparison. The total emitted energy is 
fixed, but the S5 simulation has more energetic photons 
which penetrate the ionization front further than those 
of S4. However, the difference of the gas temperature be- 
tween two simulations is negligible. The temperature of S5 
is slightly higher than S4, but the difference is less than 
0.5K at all rcdshifts. 



We estimate that our 1% model yields ar ound 5 times 
more X-rays than the fiducial model in ISantos et al 



(2008), although the source formation modeling is quite 
different and the comparison is difficult (we do not use the 
dark matter halo mass at all in computing the star forma- 
tion rate). However, comparing their plots of the average 
temperature and ionization fraction evolution with ours, 
we can deduce that their average gas kinetic temperature 
rises above the CMB temperature around ionization frac- 
tion Xbii = 10% while in our case the same event occurs at 
xmi = 15%. We find several reasons for this apparent dis- 
crepancy. First we defined neutral IGM as xhii < 0.01 
and use this to compute the average gas temperature. 
Although this is not absolutely explicit in the paper, we 
believe they use xhii < 0.5 , thereby including warmer gas 
in the average. Then, they have a more extended reion- 
ization history, which reduces the effect of the delay in 
the X-ray heating (see next section). Finally the initial 
X-ray heating is shifted to higher redshifts, when the dif- 
ference between the average neutral gas temperature and 
the CMB temperature is less. 

4.3. Brightness temperature maps 

We have run Lyman-a simulations as a further post- 
treatment to obtain the differential brightness tempera- 
ture <5Tf,. The 5Tb is determined by va rious elements, and 



it is expressed as (jMadau et all 119971 ): 



ST b « 28.1 mK x m (1 + 5) 



1 + z 
10 



To - T, 



CMB 



T s 



(7) 



where 8 is the baryon over density, T s is the spin tem- 
perature, Tcmb is the CMB temperature, and xbi is the 
neutral fraction. The contribution of the gradient of the 
proper velocity is not considered in this work. 

The spin temperature T s can be computed with: 



1 s 



T, 



CMB 



Xq/T c -\- x c T 



and 



27AioTcmb 



1 



and 



^ioTcmb 



(8) 



0) 



where P a is the number of Lyman-a scatterings per atom 
per second, A\q is the spontaneous emission coefficient of 
the 21 cm hyperfine transition, is the excitation tem- 
perature of the 21cm transition, and Cio is the deexcita- 
tion rate via collisions. Details on derivi ng these relations 
and co mputing Cio can be found, e. g., i nlFurlanetto et al.1 
(12006). The peculiar vel o city g radients (jBarkana fc Loebl . 
2005; B haradwai &: Ali . 2004) is not considered in this 
work. 

As we can see in eq. Ts is coupled to the CMB tem- 
perature Tomb by Thomson scattering of CMB photons, 
and to the kinetic temperature of the gas Tk by collisions 
and Ly-a pumping. Coupling by collisions is efficient only 
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at z > 20, or in dense clumps, so Ly-a is the key coupling 
process in the diffuse IGM. The 8Tb maps are a good way 
to see how these different elements affect the signal. In 
Fig. 2] we show several 8Tb maps of the same slice from 
different radiative transfer simulations, SI, S2, S6 and S7. 
S4 and S5, which sets the X-ray luminosity at 0.1 % of 
the UV luminosity show a trend very similar to SI and 
are not plotted. The bandwidth of the slice is 0.1 MHz for 
all maps, which corresponds to 1.9 Mpc for the maps on 
the left column (a)-(d), 1.8Mpc for (e)-(h) and 1.6 Mpc 
for (i)-(l). 

The left four maps of 8Tb in Fig. HJ (a)-(d), are plot- 
ted when the mass averaged Ly-a coupling coefficient x a 
is (x a ) = 1. This value is interesting because in this 
moderate coupling regime, fluctuations in the Ly-a lo- 
cal flux induce fluctuations in the brightness tempera- 
ture, which is not the case anymore when the coupling 
saturates. The corresponding redshifts are z — 10.50 for 
S2 and 10.13 for the others. The corresponding averaged 
ionization fractions are 0.005, 0.018, 0.005 and 0.005 for 
(a)-(d). Indeed, the averaged ionization fraction of S2 is 
higher than the others since it uses a harder spectrum. 
The ratio of the integrated energy emitted in the Lyman 
band (Ly-a < E < Ly-limit) with respect to the ionizing 
band, /3 — E^y man j 'Ei on , is three times less for S2 than 
for the others. For a given number of emitted Ly-a pho- 
tons, a harder spectrum produces a larger number of UV 
ionizing photons, therefore S2 has a higher ionization frac- 
tion when (x a ) = 1. SI shows a deeper absorption region 
around the ionized bubbles than S2, and it is also due to 
the different ratio of the number photons in the Lyman 
band and the ionizing band. In the case of SI, the ionized 
bubble is smaller than the highly Ly-a coupled region. 
Since the kinetic temperature outside of ionized bubble 
is a few kelvin, which is lower than the CMB tempera- 
ture (Tomb ~ 30 if), the neutral hydrogen has a strong 
21-cm absorption signal. On the other hand, the highly 
Ly-a coupled regions in S2 mostly resides in the ionized 
bubbles, which are bigger than in SI. S7 has almost the 
same averaged ionization fraction and ionized bubble size 
as SI, but the gas around the ionized bubbles as well as in 
the void is heated by strong X-rays. The signal is still in 
absorption because the X-ray heating has not been able 
to raise the IGM temperature above Tcmb, but the inten- 
sity is reduced. Contrary to SI and S2, the neutral gas 
around the ionized bubbles produces a weaker signal than 
in the void, because the gas around the bubbles is more 
efficiently heated by X-rays. S6 shows an intermediate be- 
havior between SI and S7. 

The four maps in the middle of Fig. (e)-(h), are 
for (x a ) = 10. The corresponding redshifts are z — 9.03 
for S2 and 8.57 for the others. The averaged ionization 
fractions are 0.043, 0.141, 0.043 and 0.040 for (e)-(h). 
These redshifts when (x a ) = 10 are interesting because 
the amplitude of fluctuation 8Tb reaches a maximum. If 
we do not consider the effect of Ly-a coupling and as- 
sume Tk ^> Tcmb, which does not allow the signal in 
absorption, the largest fluctuations would appear around 



{xmil 



.5 as n oticed by iMellema et al. ( 2006aF) and 



Lidz et al.l (|2007al ). but including the inhomogeneous Ly- 
a coupling and computing Tk self-consistently, this max- 
imum is shifted to an earlier phase of reionization. The 
ionization fraction and bubble size in S2 are still greater 
than in the other models, but the absorption intensity is 
lower than in SI. Here is why. The evolution of kinetic 
temperature in the void regions is dominated by the adia- 
batic cooling: the temperature drops as the expansion pro- 
gresses. The kinetic temperature of SI in the voids is lower 
than in S 2 by 0.5 K due to the difference in redshift, which 
explains the stronger absorption intensity in SI. The neu- 
tral gas around ionized bubbles in S7 is heated by the high 
X-rays level above Tcmb , and starts to produce the signal 
in emission. The neutral gas in the voids is also affected by 
X-rays. However, it is not sufficiently heated yet so that 
the signal turns everywhere from absorption to emission. 
Nevertheless the intensity of the signal is reduced by the 
X-ray heating, and it shows the weakest signal among the 
four maps at (x a ). 

The four maps on the right of Fig. 0] (i)-(l), are for 
(zhii) = 0.5. The corresponding redshifts are z = 7.68 for 
S2, 7.00 for SI and S6, 6.93 for S7. The averaged Ly-a cou- 
pling coefficients, x a , are 138.3, 34.9, 138.75 and 187.75 
for (i)-(l). Contrary to the above cases, the absorption in- 
tensity of SI in the void region is weaker than that of S2. 
This is due to the Ly-a heating. Ly-a heating is negligible 
during the early phase, but the amount of Ly-a heating 
accumulated between z ~ 12 and z ~ 7 can heat the 
gas in the voids by several kelvins. In order to reach 50% 
ionization, SI produces a larger number of Ly-a photons, 
which propagate beyond the ionizing front. The (x a ) of 
SI is almost 4 times greater than that of S2, and the ac- 
cumulated Ly-a heating increases the kinetic temperature 
by 3-4 K more than in S2. The intensity in absorption is 
very sensitive to the value of kinetic temperature, so this 
small amount of heating reduces the signal by up to 100 
mK. The X-ray heating in S 7 is strong enough to heat all 
the gas above Tcmb at this redshift, so we see the signal in 
emission everywhere. S6 shows intermediate features be- 
tween SI and S7, showing the weakest signal. Indeed, the 
X-ray heating in S6 increased the gas temperature in the 
neutral voids just around the Tcmb, which is the transi- 
tion phase from absorption to emission. 

4.4. Power spectrum 

Fig.[5]shows 3 dimensional power spectra of the brightness 
temperature fluctuations for the SI, S2, S6 and S7 simula- 
tions. The power spectrum can be defined as the variance 
of the amplitude of the Fourier modes of the signal for a 
given wavenumber modulus: 



P(k) = (8T b (k)8T b (-k)>. 



(10) 



We binned our modes with 8k = ^(Mpc/ft.) -1 and plot- 
ted the quantity A 2 = k 3 P(k)/2TT 2 . The power spectra 
of S4, S5 are not presented in Fig. [5] since their pat- 
terns are similar to SI. During the early phase, when 
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(x a ) = 1, the amplitude of the powerspectrum of 4 sim- 
ulations are similar. The spectra follow patterns similar 
to the power spectrum of S100 in Paper I (see Fig. 15). 
M odel S6, show s a spe ctrum similar to the fiducial model 



of Santos et al 



(120081) . The main difference in shape ap- 
pears at k > lh _1 .Mpc in our model. This is possibly 
connected to the fact that they assume a -3 dependence 
of the Lyman-a flux, while at short distances from the 
sources ( < 10 comoving Mpc), wings scatteri ng effects 



produce a dependence ( Semelin et al. . 20071) . Also no- 



ticeable is the difference between S6 and S7. S7 is depleted 
at small scale, the strong X-ray heating damping the ab- 
sorption near the sources. On very large scale, however 
the already strong heating in S7 creates temperature fluc- 
tuations which boost the S7 power spectrum. 

When (x a ) = 10, the power of both S6 and S7 de- 
creases since X-ray heating prevents a strong absorp- 
tion signal. The strong X-ray heating of S7 increases the 
gas temperature around Tcmb, and it shows the small- 
est power. However, the power of S6 falls down under 
the power of S7 when the hydrogen is 50 % ionized. At 
this redshift, the rising temperature of neutral gas reaches 
Tqmb in S6, while it is already much greater than Tomb in 
S7. This is also visible in Fig. [4] Later, when (xhu) = 0.9, 
all p ower spectra drop. Our S6 model agrees quite well 
with lSantos et al] (|2008l ) both in shape and amplitude for 
these two last stages. Indeed both the effects of Lyman- 
a coupling and X-ray heating reach a saturation in the 
determination of the brightness temperature, erasing the 
differences in our treatments. S2 has the largest power 
over all scale when (xhii) = 0.5 and (xhii) = 0.9. This 
is due to the near lack of Ly-a heating. Let us mention 
however that some sort of transition to Pop II formation 
should have occurred by then, providing some level of Ly- 
a heating. So S2 is probably not realistic during the late 
EoR. 

In brief, the 21-cm power spectra of our models vary 
in the 10 to 1000 mK 2 range, in broad agreement with 
Santos et al. (2008) who included the inhomogeneous X- 



ray and Ly-a effect on the signal in a semi-analytical 
way, with moderate discrepancies at high-redshift and 
small scale due to wing effet in the Lyman-a radiative 
trans fer. Quite logically our results differ at high-redshift 
from dMellema et al.L[2006bHZahn et alll2007tlLidz et al' 



l2007btl iMcQuinn et all 120061 ) who focused on the emission 
regime. These authors found a flattening of the spectrum 
around (xhii) = 0.5. It is interesting to notice that in the 
case of a strong X-ray heating (model S7) the spectrum is 
quite flat at all redshift (temperature fluctuation boost the 
power on large scales at high-redshift). In the future obser- 
vations, this would be a first clue of larger-than-expected 
contribution from X-ray sources. 

We now plot the evolution of the power as a function 
of redshift for 4 different k values. The evolution of the 
power spectrum with and without X-rays is very differ- 
ent. SI and S2, which do not have X-rays, show a single 
maximum on small scales (k — 1.00 h/Mpc and k = 3.15 
h/Mpc) around redshifts 8.5 and 9 which correspond to 



r — si — s2 

56 S7 <x„>=l. 


<x a >=10- 

-+++H 1 1 1 H-l-H 


-h-hH 1 1 1 I 1 1 | 













Fig. 5. Power spectrum evolution of the 6Tb from the SI 
(thin black), S2 (thick black), S6 (thin gray) and S7 (thick 
gray) simulations. S2 use top-heavy IMF whereas the oth- 
ers use Salpeter IMF. S6 has 1% and S7 has 10% of total 
luminosity for X-rays. 

the redshifts of (x a ) = 10 for each simulation. On large 
scales (k = 0.07 h/Mpc and k = 0.19 h/Mpc) the power 
spectrum shows two local maxima on large scale. The first 
peak is related to the Ly-a fluctuations. The 6Tb fluctua- 
tions are dominated by Ly-a fluctuations at high-redshift, 
but it decreases when the Ly-a coupling saturates. Then 
it rises again. This time, the fluctuations are dominated 
by the fluctuations of ionization fraction. The second peak 
appears at the redshift when (xhii) = 0.5 for each simu- 
lations. The overall amplitude of SI and S2 are similar, 
but the position of the local maximum peaks of S2 are at 
higher redshift due to the faster reionization. The key to 
the single-double peak difference is that the contribution 
of the ionizing field fluctuation to the brightness tempera- 
ture power spectrum increases during reionizatio n on large 
scale but not on small scales ( Iliev et al 1 l2006bh . 

With X-rays (models S6 and S7), the evolution fol- 
lows a different sc enario. We find a pattern similar to 
Santos et all ([20081 ). On small scales (thin and thick gray 
in Fig. [6]) , the intensity of the signal increases up to the 
maximum as the spin temperature couples to the kinetic 
temperature. Then it decreases during the absorption- 
emission transition. As the fluctuation due to the ioniza- 
tion fraction comes to dominate, the power reincreases 
slightly or remains in a plateau until it drops at the end 
of reionization. The evolution of the power on small scale 
does not show a marked minimum. The evolution of large 
scales (thin and thick gray in Fig. ^ is the most inter- 
esting: it shows three maxima. From high-redshift to low 
redshift, each peak corresponds to the period where the 
fluctuation of the Ly-a coupling, gas temperature and ion- 
ization fraction dominate. There exists a deep suppression 
between the second and the third peaks which does not 
appear without X-ray. It occurs when the X-ray heating 
raises the gas temperature of the neutral IGM around 
TcmBi which dampens the signal. The second minimum 
in S7 occurs earlier than in S6, since the stronger X-ray 
heating of S7 increase the gas temperature around Tqmb 
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at a higher redshift. We find a much narrower third peak 
in S4 (not plotted), which uses a 10 times weaker X-ray 
heating than S6. The position and amplitude of the peaks 
as well as the width depend on the intensity of X-ray heat- 
ing. The width of the third bump (6.5 < z < 7.5) is the 
largest in S7 (with 10% X-ray) and the smallest or negli- 
gible in S4 (with 0.1% X-ray). The existence/position of 
this third peak and of the second dip in the evolution of 
large scale power spectrum will be measurable by LOFAR 
and SKA observations, and it will help us constrain the 
nature of the sources during the EoR. 
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Fig. 6. Evolution of the brightness temperature power 
spectrum with redshift. k = 0.07 h/Mpc (thin black), 
k = 0.19 h/Mpc (thick black ), k = 1.00 h/Mpc (thin 
gray) and k = 3.15 h/Mpc (thick gray) 



4.5. Non-Gaussianity of the 21-cm signal 



The non-gaussianity 
studied in pr evious 



of the 
works. 



21-cm signal h as bee n 



Ciardi fc Madaul (120031) : 



Mellema et al. (l2006bl) shows the non-gaussianity of the 
21-cm signal in numerical simulati ons by computing the 
Pixel Distribution Function (PDF). lichikawa et all (|200°J) 
draws the histor y of reionization from the measurement of 
the 21-cm PDF. Hark er et al. I (|2009l) compute the skew- 
ness of the PDF and show how it could help in separating 
the cosmological signals from the foregrounds. However, 
all these works model the signal in emission only. 

We present the 21-cm PDF from our simulations in 
Fig. [71 for several representative redshifts. To obtain 
the 21-cm PDF, we sample 5Tb within a 1 h~ 1 Mpc 
resolution, an acceptable value for SKA. The 21-cm 
PDF from our simulations is highly no n-Gaussian as ex- 
pected , but is also qu it e diffe re nt from Ciardi fc Madaul 
(120031): iMellema et all (j2006bl) ; lHarker et al.1 (|2009f) and 
Ichikawa et al.1 ( 2009 ). Our distributions extend to nega- 
tive differential brightness temperature with a variety of 
shapes depending on the redshift. 

The panel (a) of Fig.[jjis the 21-cm PDF at the begin- 
ning of reionizati on, z = 14.05. It is a t the beginning of 
reionization that Ichikawa et al. (2009) find a 21-cm PDF 



closest to a Gaussian, but this is also when their model 
is the less relevant. In our case, all signals are found in 
absorption and their distribution is peaked around OmK, 
completely non-Gaussian. 

We show the PDFs at z = 10.64, still during the early 
EoR, in panel (b) of Fig. [71 The position of peaks are 
shifted around -lOOmK ~ -50mK which means that the 
spin temperature of the particles is decoupled from Tomb 
by Ly-a photons. The PDF is much close to a gaussian dis- 
tribution, extending to positive values (the smooth curves 
are the best gaussian fit to the PDFs). However, all of 
them are left skewed (toward negative temperature). The 
reason for this negative skewness is the same as the rea - 



son for the po s itive s k ewness in Mellema et al. (2006b) 



Ichikawa et ail (|2009l) : barker et alj ([20091 ): it is due to 



the signal from high density regions seen in absorption for 
us and in emission for them. 

Panel (c) in Fig. [7J shows the PDFs when z = 8.48. 
Here we find a bimodal distribution with a plateau in the 
absorption region between 100 mK and mK, for models 
SI, S2 and S4. The left peaks of the PDFs moves also 
toward higher 5Tf> with increasing heating efficiency. In 
the case of S7, the left peak merges with the right one, 
and the form is very similar to a gaussian. 

The panel (d) of Fig. [7J is plotted when the ionization 
fraction is 50%. The width of the PDF of SI and S4 is re- 
duced because Ly-a heating is well advanced. We find sig- 
nals in emission in S6 and S7, since X-rays heat the gas in 
neutral re gions above Tqmb- Ind eed, these PDF forms are 
similar to lichikawa et all (|2009l ). The PDFs of S6 and S7 
could be fitted by the Dirac-exponen tial-Gaussian distri- 
bution used bv lichikawa et all ([20091 ) . S2 has broad PDF 
still, since the Ly-a heating is 4 times lower than in the 
others and it retains the signal in strong absorption. 

It is interesting to note that the PDF always shows 
a spike around ~ mK. During the beginning of reion- 
ization, it is due to the large amount of neutral hydrogen 
whose spin temperature is still well coupled to the CMB 
temperature. As the reionization proceed, the spin tem- 
perature is decoupled from the CMB and the number of 
pixels at ~ mK decreases, but the peak grows again 
with the increasing contribution from completely ionized 
regions. This feature is interesting since interferometers 
such as LOFAR or SKA only measure fluctuations in the 
signal and do not directly provide a zero point. 



Ichik awa et al.1 (|2009|) extract informations about the 
averaged ionization fraction from the 21-cm PDF. As 
could be expected our PDFs converge with their result 
when the contribution from X-ray sources is sufficient and 
the EoR somewhat advanced. What we find is that a clear 
tracking of the nature of the ionizing source remains on 
the PDF when the absorption phase is modeled. A strong 
X-ray contribution produce unimodal PDFs while a weak 
X-ray contribution yield a bimodal PDF. 
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The evolution of the skewness is presented in Fig. \E\ 
The skewness 7 is defined as 



7 



ST b f 



Hi(6T^-5T b ) 2 ]V 2 



(11) 



LJV 



where N is total number of pixels in 5Tt data cube, ST^ is 
STi, in i th pixel, and 8Tb is the average on the data cube. 

At the beginning the skewness is highly negative for 
all simulations as we can expect from the panel (a) and 
(b) of Fig. [71 Then, in all models, the skewness rises to a 
local positive maximum when the average ionization frac- 
tion is a few percents. It is interesting to notice that the 
skewness of all simulations close to zero again when the 
neutral fraction is about 0.3. While this behavior could be 
used to provide a milestone of reionization, its robustness 
should first be checked. We can notice that two of the 
three models presented in lHarker et al. show the 

same behavior. Those are however the two less detailed 
models. 

Another interesting feature in Fig. [5] is that the skew- 
ness of S7 has two local maxima while others do not. 
Again, this could be used as a clue to a large contribu- 
tion from X-ray sources. 





1 1 1 1 1 1 1 1 1 




1 ■' 
/ ' 

/// 










/ 1 1 
' 1 ■' 
/ / 










*• ■ * 










/ 

/ 

/ 

/ 










/ 








— SI 










S2 










-- S4 
-■■ S6 
■-■ S7 








I , , , 1 I 1 I 1 1 1 , , , 1 , , I I 1 1 1 1 1 1 , , , , 1 , I I 




1 1 t t t t I 1 i l 1 





0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

Ionization fraction <x,,,, > 



Fig. 

tion. 



8. Evolution of the skewness with the ionization frac- 



5. Conclusions 

We modeled the 21-cm signal during the EoR using nu- 
merical simulations putting the emphasis on how various 
types of sources can affect the signal. The numerica l meth - 
ods used in this work are similar to Baek et al. (|2009h . 
The N-body and hydrodynamical simulations have been 
run with GADGET2 and post-processed with UV contin- 
uum radiative transfer further processed with Ly-a trans- 
fer using LICORICE, allowing us to model the signal in 
absorption. The main difference from the previous work is 



a more elaborated source model, including X-ray radiative 
transfer and He chemistry. 

We have run 7 simulations to investigate the effect of 
different IMFs, helium, different spectral indexes and the 
different luminosities of X-rays sources. The reference sim- 
ulation in this work, SI, using only hydrogen and stellar 
type sources, reached the end of reionization at z » 6.5 
and showed a strong absorption signal until the end of 
reionization. Our top heavy IMF (model S2) produces 
~ 2.6 times more ionizing photons than the Salpeter IMF. 
S2 reached the end of reionization earlier than the others 
by Az w 1. In addition the different SED changes the ratio 
of Ly-a and to ionizing UV photon numbers, and it slows 
down the saturation of the Ly-a coupling and the heating 
by Lyman-a in the top heavy IMF case. This modifies the 
statistical properties of 21-cm signal. 

The simulation with helium, S3 also has a slightly ear- 
lier reionization than the others since the number of emit- 
ted photons per baryon is higher. Except for the slightly 
lower kinetic temperature in the bulk of ionized regions 
due to the higher ionization potential than for hydrogen, 
the properties of the 21-cm signal from S3 is similar to SI. 
We chose QSO type sources with a power-law spectrum 
as X-ray sources in model S4 to S7. The spectral index 
a has large observational uncertainty, so we used two dif- 
ferent spectral indexes. S4 and S5 have 0.1% of the total 
luminosity in the X-ray band. S5 uses a = 0.6, while other 
simulations with X-rays use a = 1.6. S5 showed very little 
difference on the gas temperature with respect to S4. S4, 
S6 and S7 have different luminosities in the X-ray band, 
keeping the same values for the other simulation param- 
eters. Using a stronger X-ray luminosity indeed increased 
the gas temperature in the neutral hydrogen. Accordingly 
the 21-cm signal and its power spectra are modified. We 
found an increase of a few kelvin for the neutral gas tem- 
perature in our fiducial model, S4, in which X-rays account 
for 0.1% of the total emitted energy. The 21-cm signal in 
S4 was similar to SI, showing the maximum intensity in 
absorption, ~ 200 mK, at z s» 9. Stronger X-ray levels 
increase the gas temperature and reduce the intensity. We 
found that in S6 and S7, which uses 1% and 10% of the 
total luminosity for X-rays, the absolute maximum inten- 
sity in absorption decreases to ~ 130 mK and ~ 80 mK. 
The 21-cm power spectrum of our work is greater by two 
or three orders of m agnitude than in works focusing on 
the emission regime dMellema et all l2006bl IZahn et al ' 



2007t iLidz et all l2007ri iMcQuinn et all 120061) . However. 



the results are in broad agreement with the work of 
ISantos et al. (2008), who modeled absorption using semi- 
analytical methods for X-ray and Lyman-a transfer. We 
noticed that the 21-cm fluctuation is dominated by Ly-a 
fluctuations during the early phase, X-rays later (or the 
gas temperature), and the ionization fraction at the end. 
This is visible on the evolution of the 21-cm power spec- 
trum with redshift. The 21-cm PDF of our work was dif- 
ferent from other work, since we do not assume that the 
spin temperature T s ^> Tcmb- 



S. Baek, B. Semelin, P. Di Matteo, Y. Revaz, F. Combes: Reionization by UV or X-ray sources 



13 



The first most important conclusion from our work is 
that even including a higher than generally expected level 
of X-ray, the absorption phase of the 21-cm survives. Its 
intensity and duration are reduced, but the signal is still 
stronger than in the emission regime. Heating the IGM 
with X-rays takes time! 

The second important result is that we found three di- 
agnostics which could be used in the analysis of future ob- 
servations to constrain the nature of the sources of reion- 
ization. (i) The first and maybe the most robust is the evo- 
lution with redshift of large scale modes ( k ~ 0.1 h/Mpc) 
of the powerspectrum. If reionization is overwhelmingly 
powered by stars, this evolution should have one local 
minimum (two local maxima) . However, if the energy 
contribution of QSO is greater than ~ 1%, a second lo- 
cal minimum (third maximum) appears. The higher the 
X-ray level, the broader the third peak, (it) The second 
simple diagnostic is the bimodal aspect of the PDF which 
disappears when the X-ray level rises above 1% of the to- 
tal ionizing luminosity. (Hi) Last is the redshift evolution 
of the skewness of the 21-cm signal PDF. While all other 
models show a single local maximum at a few percent 
reionization, a very high level of X-rays (> 10% of the to- 
tal ionizing luminosity) produces a second local maximum 
appear around 50% reionization. 

Modeling the sources in the simulation is complex. It 
involves taking the formation history, IMF, SED, life time, 
and more into account. Although detailed models are de- 
sirable for the credibility of the results, we believe that the 
effect in the 21-cm signal can be bundled in 3 quantities. 
The first is the efficiency: how many photons are produced 
by atoms locked into a star. This parameter must be cal- 
ibrated to fit observational constraint: end of reionization 
between redshift 6 and 7, and Thomson scattering opti- 
cal depth in agreement with CMB experiment. The two 
other quantities which contain most of the information are 
two box-averaged ratios: the energy emitted in the Lyman 
band to the energy emitted in the ionizing band ratio and 
the same ionizing UV to X-ray ratio. In this work we ex- 
plored values of 0.32 (model S2) and 0.75 (all other mod- 
els) for the former and 0.001, 0.01 and 0.1 for the latter. 
Once additional physics is included in the simulation and 
using a higher resolution to account for all the sources, it 
will be interesting to explore the value of these quantities 
systematically. 

We mentioned in the introduction that the minimum 
boxsize for reliable predictions of the signal is 100 /i _1 Mpc. 
It is important to realize that t his value (confirm ed by 
emission regime simulations, e. g., Iliev et al. 2006bl ) is es- 
timated based on the clustering properties of the sources 
and applies to the topology of the ionization field. It 
may be underestimated when we study the early absorp- 
tion regime, when only the highest density peaks con- 
tain sources. Their distribution is the most sensitive to 
possible non-gaussianities in the matter power spectrum. 
Moreover, they are distant from each other and, conse- 
quently, produce large scale fluctuations in the local flux 



of Lyman-a and X-ray photons. We intend to extend our 
investigation to larger box sizes in a future work. 

A few final words on additional physics not included 
in our model. Shock heating from the cosmological struc- 
ture formation is ignored, but it could have the poten- 
tial to affect the 21-cm signal by increasing the gas tem- 
perature above the CMB temperature. However, it is not 
sure whether shocks are strong enough in the filaments of 
the neutral regions to affect the 21-cm signal. Mini-halos 
(~ 10 4 — 10 s Mq) form very early during the EoR and are 
dense and warm enough from shock he ating during viri- 
alization to emit the 21 cm signal, but iFurlanetto fc Oh 
(200G) find that the contribution of mini-halos will not 
dominate, because of the limited resolution of the instru- 
mentation. However, shock heating is worth investigat- 
ing with coupled radiative hydrodynamic simulations with 
higher mass resolution. Also worth investigating is the ef- 
fect of including higher Lyman lines in the radiative trans- 
fer. 
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Appendix A: Implemented method for UV and 
X-ray 

In this section, we explain the numerical methods that 
we use for radiative transfer with helium and X-ray. A de- 
sc ription of the main methods used in LICORICE appears 



Baekl ([2009) 0. 



A.l. Hydrogen and helium ionization 

The radiation field is discretized into monochromatic pho- 
ton packets which are emitted with random direction 
and frequency (with appropriate distribution) from the 
point sources. Photon packets propagate through radiative 
transfer cells and deposit photons and energy depending 
on the absorption probability of each cell. We need to com- 
pute absorption probabilities for each absorbers H°, He 
and He+0 

The probabilities of the photon being absorbed by each 
elements are given by 



T H 



T H° + T Hc° + r Ho+ 

. ^HeO 

T R° + T Ho° + T Ho+ 
r Hc+ 



(l-c- T ), 



(A.l) 

(A.2) 
(A.3) 



7h° + T"Ho° + r Hc+ 

where rjjo, r Hc o and rjj+ are the optical depths for H°, 
He and He + in a given cell, and r = t h o + t Hc o + r H + . 
Therefore, the total absorption probability, Vtotai, for a 
photon packet arriving in a cell of optical depth r is given 

by 

Vtotaiij) = 7> H « + V Hc o + 7> Hc + = 1 - e" T . (A.4) 
For each cell, the optical depth r is 



t = r H o + r Hc n - 



T Hc+ 

- O-He (v)n Hc o + (T Hc +(^)"-Ho+] h ( A -5) 



where a a is the photoionization cross-section for absorber 
A £ {H°, He , He + } , n l A , is number density in the cell and 
I is the path l ength in a cell . We u se photoionization cross 
section fits in Verner et al. I (Il996h . 

If _/V 7 is the number of photons in a photon packet 
arriving in a cell, then the number of photons absorbed in 
a cell, Na, is : 



N A = N~V, 



7 ' total 



AT 7 (l 



(A.6) 



Available at http://aramis.obspm.fr/~baek/these.pdf 

7 H° denotes the neutral hydrogen and He denotes the neu- 
tral helium. 
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The contribution of photoionization and photoheating to 
the evolution of the ionization fraction and the tempera- 
ture within a time step At is: 

(A.7) 



Ax H + = — r H °At = —!- — 
n H iV H 



AT = 



A* Ho p *, N y P Hc o 
Ax He + = 1 He At = — , 

At "Hep a, ^V Ho+ 

A.t Ho+ + = 1 Hc + At = — 

riH A Ho 



-| - ^k B TAn + n~ ( Vna{hv - hv thjH o) 



(A.8) 
(A.9) 



+ n n V ne a(hv - hv thMc a) 



(A.10) 



where TVh and Nn e is the total number of hydrogen and 
helium in a cell, hv t h,A are the ionization potential of the 
recombined atom, n 7 is the iV 7 in a unit volume. V a are 
the continuous photoionization rates used to actually in- 
tegrate the evolution of th ionization fractions. Indeed re- 
combination, collisional ionization and radiative cooling 
are treated as continuous process, with integration time 
step At* much less than At using the following coupled 
equations 

"H dX , H+ = 7 H o(r)n H on e - a H +(r)n H +n e + r H nn H o , 
= lHc°(T)n Hc on e - 7 Ho +(T)n Hc +n e 



dt 
dx He + 



dt 



-a Ho +(T)n Ho +n e (A.ll) 
+a Rc+ +(T)n Hc ++n e + r Hc on Ho o , 



n He 7, = 7He+U J^Hc+ne ~ «Hc++ U )%e++"e 



dt 



_ rHo+ n Hc+: 



where a/ and ^a indicate the recombination and colli- 
sional ionization coe fficients and I g (H + , He + , He ++ | . 
We use the values in Hui fc Gnedinl (|l997h for the recom- 
bination coefficient and recombination cooling. For colli- 
sional ionization co efficients an d other radiative cooling 
we use the values in ICenl (|1992| ) . 

The temperature is computed from the energy con- 
servation equation. E is the internal energy of the gas, 
E = ^nk B T. 



T-L and A are the heating and cooling function which ac- 
count for the energy gained and lost in a unit volume per 
unit time. 

Each photon packet keeps propagating until it exits 
the simulation box (if we use free boundary condition) or 
until the remaining photon content is much less than the 
initial photon content, N 1 < 10~p N^ nlUal . We typically 
adopt p = 4 for the UV continuum. 



Appendix B: Validation tests 

B.l. Radiative Transfer Comparison Test 

A cosmological r adiative transfer co de comparison project 
was performed in llliev et al.l ( 2006a ) trying to understand 
which algorithms (including various flavors of ray-tracing 
and moment schemes) are suitable for a given non-trivial 
problem as well as to validate each code by comparing the 
results with other codes. Five tests are run for radiative 
transfer in a static density field. We reproduced several of 
tests in llliev et al.l (|2006al ) with LICORICE. LICORICE 
shows good agreements with other code s: resu l ts and com- 
parisons with other codes are shown in iBaekl (|2009l). 
Three additional tests are presented in llliev et al 



(2009) for the coupled gas dynamical and radiative trans- 
fer evolution. LICORICE directly participated in this sec- 
ond project and shows good agreements with the other 
codes. 

B.2. Comparison with CLOUDY I : helium 

An analytic solution to the radiative transfer evolution 
of an homogeneous medium around a single source exists 
only in the isothermal case with only hydrogen. In or- 
der to validate helium ionization and spectr um hardening, 
we rep roduced the Stromgren sphere test in iMaselli et al. 
(|2003l ). and compared our results with the 1-D radiative 



transfer code CLOUD^O (C08 version of the code). 

A point source, emitting as a black body at T — 60000 
K with ionizing luminosity L = 10 38 ergs _1 , is located at 
the center of the simulation box in a homogeneous den- 
sity field with n — lcm~ 3 composed of hydrogen (90% by 
number) and helium (10% by number). The gas is initially 
completely neutral at a temperature of T = 10 2 K in the 
entire simulation box cube of L\, ox — 128 pc. The com- 
parison is performed at a time t s = 6 x 10 5 yr « 5tf ec 
(where t^ ec is the characteristic time scale for hydrogen 
recombination) . 

The script used for CLOUDY is the following: 

blackbody, T=60,000K 
luminosity 38 

radius 0.01 60 linear parsecs 
hden -0.0457 
abundances all -15 
sphere static 

element abundance helium -0.9542 
punch element hydrogen "hydrogen_caseA.dat" 
punch element helium "helium_caseA.dat" 
punch temperature "temperature_caseA.dat" 
punch continuum units eV "spectrum_60pc.dat" 



In Fia lB.ll and Fig lB.21 the comparison between 
CLOUDY (solid lines) and LICORICE (circles) is shown. 
The value of the different physical quantities is plotted 
as a function of the distance from the source, expressed 
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in cell units, Aa; = lpc. The points represent spherically 
averaged LICORICE outputs. LICORICE shows a good 
agreement with CLOUDY. The positions of various ioniz- 
ing fronts agree within 1%. To obtain this good agreement 
we consider the effect of secondary ionization and heating 
by fast electrons as CL OUDY does. The secondary ion iza- 



tion and heating fits in lShull fc van Steenberd (|l985h are ^ 



only accurate for the prim ary photoelectrons higher than 
100 eV , so we follow fits in iFurlanetto fe Johnson Stoever 
(|2010h which is valid for ones less than 100 eV. 



The temperature profile shows also a good agreement, 
except for the warm tail extending beyond the ionizing 
front for CLOUDY. This warm tail may originate from 
heat conduction which LICORICE does not implement. 




p 10000 
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Fig. B.2. Comparison of temperature distribution be- 
tween LICORICE (circles) and CLOUDY (solid lines), as 
a function of distance from the point source in cell units 
(1 pc). 
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14 pc CLOUDY 
40 pc CLOUDY 
51 pc CLOUDY 
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Fig. B.l. Comparison of different ionization fractions be- 
tween LICORICE (circles) and CLOUDY (solid lines), as 
a function of distance from the point source in cell units 
(1 pc). 



Treating properly the spectrum hardening is an impor- 
tant issue in the radiative transfer with helium. Since each 
absorber (H°, He and He + ) has a different ionizing poten- 
tial, the emitted spectrum will be strongly depleted just 
above the different ionizing frequencies once it hits a suf- 
ficient amount of absorbers. We plot in Fig |B.3l the lumi- 
nosity per unit surface at different distances and compare 
between CLOUDY (solid lines) and LICORICE (circles). 
Our results agree well with the ones of CLOUDY. The 
noise at higher frequencies (> 5 Ryd) is due to the nature 
of the Monte Carlo method: the tails of the distributions 
are poorly sampled. At 4 pc, all species are already ionized 
so the medium is transparent and the spectrum is close to 
the initial one. At 14 pc, only photons with energy higher 
than 4 Ryd are absorbed to double ionize He. Above 40 
pc, we can observe also the depletion at ~ 1.8 Ryd, which 
corresponds to the single ionization of helium. 
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Fig. B.3. The luminosity per unit surface at different dis- 
tances. Comparison between CLOUDY (solid lines) and 
LICORICE (circles). 

B.3. Comparison with CLOUDY II : X-ray 

We have also performed a comparison test with CLOUDY 
to validate X-ray radiative transfer. The geometrical ini- 
tial set up is the same as in the test above. Now a point 
source in the center emits X-ray photons with energy from 
100 eV to 2 keV, and a power-law spectrum with a = —1.6 
in EqJS] The total luminosity is the same as in the helium 
test, L = 10 38 ergs _1 . The gas is composed of pure hy- 
drogen, and completely neutral at the beginning. The re- 
combination time is proportional to the inverse of electron 
density. But instead of a sharp ionization front we have a 
gradual decrease in the ionization fraction with increasing 
radius. At the boundary of the box, the equilibrium ion- 
ization fraction is only ~ 5% so the recombination time is 
~ 20 times large than in the inner region. Consequently 
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the integration with LICORICE was extended to t = 2500 
Myr. 

The script used for CLOUDY is the following: 

interpolate (1. -10.) (7.34 -10.) (7.35 0.) 
(147.35 -2.08) (147.4 -10.) (200. -10.) 
luminosity 38 

radius 0.1 100 linear parsecs 
hden -0.0457 
abundances all -15 
sphere static 

punch element hydrogen "hydrogen_Xray.dat" 
punch temperature "temperature_Xray.dat" 



In Fig |B.4l and Fig |B.5[ the comparison between 
CLOUDY (solid lines) and LICORICE (circles) is shown. 
The value of the different physical quantities is plotted as a 
function of the distance from the source, expressed in cell 
units, Ax — lpc. The points represent spherically aver- 
aged LICORICE outputs. LICORICE shows a very good 
agreement with CLOUDY. The (very) small difference in 
the ionization fraction at large radii is due to the fact that 
we stopped the integration at t — 2500 Myr, which is not 
much more than one local recombination time: the equi- 
librium is not fully established. We could integrate for a 
longer time, but then the problem would only be shifted 
to larger radii and smaller ionization fractions. 




Fig. B.4. Comparison of different ionization fractions be- 
tween LICORICE (circles) and CLOUDY (solid lines), as 
a function of distance from the point source in cell units 
(1 pc). 
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Fig. B.5. Comparison of temperature distribution be- 
tween LICORICE (circles) and CLOUDY (solid lines), as 
a function of distance from the point source in cell units 
(1 pc). 

The properties of a star of mass M are described by 
the bolometric luminosity L(M), the effective blackbody 
temperature T c g(M), and its life time tuf e (M). We can 
write the emissivity of a single star e„ (erg.s _1 .Hz _1 ): 



e„(M) = L(M) 



B v {T cS {M)) 
J °° B v (T eS (M))dv 



where B v is the Planck law. Then we can write the time 
dependent emissivity of the complete population e to t as: 



etot(M) 



jaM)dM 



H(tm{M) - t)e v {M)t;(M)dM 



were H is the Heavy side step function and the integrals 
over the masses are bounded by the mass interval in which 
we want to apply the IMF. With this quantity, we define 
a characteristic life-time in a spectral band [v\ , 1*2] : 



= 2 



■ C /o°°etot(y, t)tdtdv 
■C Jo°° e t°t(^ *) dtdv 



If the total emissivity was a simple decaying exponential, 
this definition would give two characteristic decay times. 
We choose it as a typical time during which most of the 
energy has been emitted. From this we can define the typ- 
ical constant luminosity emitted in the band as: 



Appendix C: Luminosity, SED, and typical life 
time of the sources 

Let us consider a source particle containing a population 
of stars. The number of stars dN in a mass interval dM 
is described by the IMF function £(M): 

dN = £(M)dM 



L„ 



1 



'1/1,1/2 J Vx JO 



e tot {v,t) dtdv 



The value of t Vi>V2 and L Vl>V2 are computed in table 3 for 
different IMF in the Lyman and ionizing bands. Finally, 
we can compute the time-averaged total emissivity: 



etotO)oc / e to t(M)cftoc t me (M)e„(M)£(M)dM 
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This gives the constant SED to use during a characteristic 
life-time, and the normalization is given by the relation: 

\ e tot (v)dv = L VltV2 
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Fig. 4. Differential brightness temperature maps for different simulations. The thickness of the slice is « 2 Mpc. The 
maps in the left column are when (x a ) = 1, the middles when (x a ) = 10 and the right when (xhii) = 0.5. The black 
contour separates absorption and emission region. (1) shows no absorption region. 




Fig. 7. The evolution of the 21-cm PDF. The redshifts are 14.05, 10.64 and 8.48 for (a), (b), and (c). The PDFs of 
panel (d) is chosen so that the ionization fraction is 0.5. The red curves on (b) are gaussian fits with the mean and 
variance of the PDFs. 



